%  Copyright (C) 2002 Regents of the University of Michigan, portions used with permission 
%  For more information, see http://csem.engin.umich.edu/tools/swmf
%\documentclass{article}
%\newcommand{\BATSRUS}{BATS-R-US}
%\newcommand{\bx}{\mathbf x}
%\newcommand{\bF}{\mathbf F}
%\newcommand{\bv}{\mathbf v}
%\newcommand{\bB}{\mathbf B}
%\newcommand{\divB}{\nabla\cdot\bB}
%
%\begin{document}

\section{Time stepping \label{section:time_stepping}}

\BATSRUS\ provides several options for the temporal discretization of
the system of MHD equations (\ref{eq:continuity_conserve})--(\ref{eq:energy_conserve}), 
which can be written in the compact form
\begin{equation}
\frac{\partial W}{\partial t} + \nabla\cdot \bF = S
\label{q-compact}
\end{equation}
where $W$, $F$, and $S$ represent the set of conservative variables, 
the corresponding fluxes and source terms, respectively. 
First of all, one can choose between time accurate mode and local time 
stepping.
The latter mode should be used to find a steady state solution, since
it converges much (10 to 15 times) faster than the time accurate mode.
Both modes can be combined with a 1-stage and a 2-stage scheme, although
we strongly suggest the use of the 1-stage scheme for local time stepping.
In addition to the above choices, some of the source terms can be treated
explicitly as well as implicitly.

\subsection{Local time stepping \label{section:local}}

Local time stepping means that each computational cell uses a time step
which is based on the local numerical stability criterion. Having different
time steps for each cell corresponds to a discretization of 
\begin{equation}
\frac{1}{a}\frac{\partial W}{\partial t} + \nabla\cdot \bF = S \label{q-local}
\end{equation}
where $a$ is the acceleration factor, which is a function of $W$. 
Clearly, if a steady state is found for the system of equations 
(\ref{q-local}), the first term must be zero irrespective of the value of 
the acceleration factor $a$. {\bf Therefore a steady state solution of
the accelerated system of equations
(\ref{q-local}) is a steady state solution of the original equations
(\ref{q-compact}) as well.} 

The same holds on the discrete level if local time stepping is combined
with a one-stage time stepping scheme. An explicit 1-stage (forward Euler) 
time discretization can be written as
\begin{equation}
W^{n+1}=W^n + \Delta t\left(-\nabla\cdot\bF^n + S^n\right) \label{q-one-stage}
\end{equation}
where $\Delta t$ is the {\bf local time step}. 
If a numerical steady state is found then $W^{n+1}=W^n$, and thus
the divergence of the numerical flux has to balance the source terms
exactly, again, irrespective of the value of the local time step.

The numerical value of the time step is determined by the 
Courant-Friedrich-Lewy (CFL) stability condition
\begin{equation}
\Delta t = C \left(\frac{c_x+|u_x|}{\Delta x}
                  +\frac{c_y+|u_y|}{\Delta y}
                  +\frac{c_z+|u_z|}{\Delta z}\right)^{-1} \label{q-local-dt}
\end{equation}
where $c_x$, $c_y$ and $c_z$ are the fast magnetosonic speeds in the
$x$, $y$ and $z$ directions, and 
\begin{equation}
   C<1
\label{eq:cfl_limit}
\end{equation}
is the {\bf CFL number}. 
Stiff source terms might pose additional stability constraints, 
although this is not typical for the applications of \BATSRUS.

\subsection{Time accurate mode \label{section:time_accurate}}

For time accurate mode the code must use the same time step in every
computational cell.
Numerical stability requires that this  {\bf global time step}
is determined as the minimum of the local time steps given by
(\ref{q-local-dt}). 
The global time step can become very small if $(c+|u|)/\Delta x$ is
very large in some cells. The time step is typically
determined by the smallest cells. 

For time accurate mode the forward Euler scheme (\ref{q-one-stage})
gives first order accuracy with respect to $\Delta t$.
The temporal accuracy can be made second order
with the two-stage time stepping scheme. The first stage uses $\Delta t/2$
as the time step, and the second stage calculates fluxes and sources
based on the first stage results:
\begin{eqnarray}
W^{n+1/2}&=&W^n + \frac{\Delta t}{2}
               \left(-\nabla\cdot\bF^n + S^n\right)  \nonumber \\
W^{n+1}&=&W^n + \Delta t
               \left(-\nabla\cdot\bF^{n+1/2} + S^{n+1/2}\right)  
\label{q-two-stage}
\end{eqnarray}
It is easy to see that the two-stage scheme is twice as expensive as 
the one stage scheme (\ref{q-one-stage}). 
It depends on the application whether it is worthwhile to 
use the two-stage scheme. If the time steps
are very small due to the stability requirement, 
the first order scheme may be accurate enough, or
in other words, the truncation errors of the spatial discretization
may be much larger than the temporal errors due to the one-stage scheme.
In such a case it makes sense to use the less expensive one-stage scheme.
For other problems the second order accuracy of the two-stage scheme may 
become important.

\subsection{Point-Implicit Source Terms \label{section:point_implicit}}

Stiff source terms can be handled with the point implicit scheme,
so they do not impose stability constraint on the time step.
The starting point is the semi-implicit discretization
\begin{equation}
W^{n+1}=W^n + \Delta t\left(-\nabla\cdot\bF^n + S^{n+1}\right) 
\label{q-point-implicit}
\end{equation}
Note that the source is evaluated at time level $n+1$ in contrast
with the explicit discretization (\ref{q-one-stage}). This is 
a non-linear equation for $W^{n+1}$. It is sufficient to solve a
linearized form using a Taylor expansion of $S$ around $S^n$, 
which leads to the linear equation
\begin{equation}
\left(\frac{I}{\Delta t} - \frac{\partial S}{\partial W}\right)
\left(W^{n+1}-W^n\right) = -\nabla\cdot\bF^n + S^n
\label{q-point-linear}
\end{equation}
for the unknown $(W^{n+1}-W^n)$, where $I$ is the identity matrix,
and $\partial S/\partial W$ is the Jacobian matrix at time level $n$. 
As long as $S$ depends on the local value of $W$ only, the linear
equations can be solved for each cell separately, which explains
why this is called the {\bf point-implicit scheme}. 
Even if $S$ contains some non-local terms, 
one can evaluate the non-local part at time
level $n$ and take the partial derivatives with respect to the local part 
of the source term only. For example
\begin{equation}
\frac{\partial (\bB\cdot\bv \divB)}{\partial \bB} \approx \bv\divB
\end{equation}
can be used. A partially implicit discretization of the source term
is still more stable than an explicit discretization.

The point implicit scheme has been written for some specific 
source terms, and in most applications it is not needed and it
should not be used. 


%\end{document}
